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The microscopic structure of ultra-thin oxide barriers often plays a major role in modern nano- 
electronic devices. In the case of superconducting electronic circuits, their operation depends on the 
electrical non-linearity provided by one or more such oxide layers in the form of ultra-thin tunnel 
barriers (also known as Josephson junctions). Currently available fabrication techniques manufac¬ 
ture an amorphous oxide barrier, which is attributed as a major noise source within the device. The 
nature of this noise is currently an open question and requires both experimental and theoretical 
investigation. Here, we present a methodology for constructing atomic scale computational models 
of Josephson junctions using a combination of molecular mechanics, empirical and ab initio meth¬ 
ods. These junctions consist of ultra-thin amorphous aluminium-oxide layers sandwiched between 
crystalline aluminium. The stability and structure of these barriers as a function of density and 
stoichiometry are investigated, which we compare to experimentally observed parameters. 


I. INTRODUCTION 

Josephson junctions: ultra-thin insulating layers sand¬ 
wiched between layers of superconducting metal, are 
the fundamental building blocks of the next generation 
of quantum electronics. Examples include supercon¬ 
ducting quantum-bitiPM®, low power Rapid-Single-Flux- 
Quantum circuit^®, Superconducting Quantum Interfer¬ 
ence Device^ and non-linear elements for single quanta 
microwave electronic^®!]. i n a ]j these cases, Joseph¬ 
son junctions provide the non-linear element that allows 
quantum effects to manifest in the voltage, current or 
magnetic flux signatures of these circuits. Recent work 
on superconducting qubitd 13 14 has shown that a key lim¬ 
iting factor in quantum electronics is the existence of loss 
mechanisms, which can be traced to material defects in 
the oxide coating (and protecting) the metallic circuits, 
as well as the oxide which forms the Josephson junction 
tunnel barrier. Recent experimental probes of so-called 
‘strongly coupled’ defects^ 17 have shown that they can 
be individually addressed and manipulated, and mostly 
likely reside within the junction 1 ®. It is therefore funda¬ 
mentally important to identify and ideally remove these 
defects as a source of loss and imperfection in quantum 
circuits. 

From a quantum simulation point of view, this is not 
a trivial problem. The oxide barrier of a junction is 
amorphous, so crystalline symmetries cannot be used to 
reduce the state space. Added to this is the far more 
fundamental issue that we currently do not understand 
what forms the defects of interest. Various microscopic 
models exist, includi ng hy drogenic dangling bond^ffiH, 
charged surface states- 3 24 and delocalisation of the oxy¬ 
gen atoms themselve d 25 l 26 l As well as atomistic models, 
a range of effective defect state models also exist such 
as phonon dressing of electronic stated 7 ', metal-insulator 
gap statesPl and Andreev bound state modelp®. One 
possible way of distinguishing between these options is 
to develop complete atomistic models of the Josephson 
junction and study the configuration of the amorphous 
layer. Forming such atomistic models using molecular 


mechanics and ab initio methods is the focus of this pa¬ 
per. 


II. THE JOSEPHSON JUNCTION FORMATION 
PROCESS 

Josephson junctions may be constructed from any 
superconducting material with any insulating or non¬ 
superconducting metal barrier to invoke a weak link cou¬ 
pling. A popular material choice involves the use of alu¬ 
minium as the superconducting material, and an amor¬ 
phous oxide layer as an insulating barrier. 

Shadow evaporation is a common technique used to 
fabricate a system such as this, where two metallic layers 
are deposited from different angles with an intervening 
oxidation step. This is usually performed using a Dolan 
bridge, which obscures part of the substrate during each 
metal deposition stepl 23; . It has more recently been shown 
that junction fabrication can be performed without the 
requirement of this bridgP®. Regardless of the process 
chosen, the oxidation of the aluminium does not result in 
a set of crystalline monolayers, but a non-uniform amor¬ 
phous layer varying in stoichiometrjl 31 ^] densitj® an d 

thickness^EUD (nominally ~ 2 mn). Although epitax¬ 
ial growth of aluminium-oxide barriers has been demon¬ 
strated^, this technique is not yet mainstream as it 
is considerably more difficult than conventional shadow 
mask evaporation. It is therefore the amorphous oxide 
formation which needs to be investigated predominantly, 
in order to obtain results from simulation which are ap¬ 
plicable to future fabrication work. 

Simulating oxide layer growth is in general a diffi¬ 
cult problem as the time scale of the oxide growth (~ 
minutes) is many orders of magnitude greater than typ¬ 
ically achievable molecular dynamics timescales (ps-ns). 
One standard approach is to perform the simulation at 
elevated temperatures and gas pressures (> 1 atm^P®®]. 
This accelerates the oxidation process, making the com¬ 
putation feasible on current high performance comput¬ 
ing infrastructure. However, it also removes the Simula- 
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tion from the reality of experimental junction formation, 
where pressures range between 10 -9 and 10 -3 atm®®* 
It remains to be seen whether any fundamental physics 
is neglected by adopting this approximation. 

An alternative approac h is to form an amorphous layer 
via direct melt and quenclP®^. This method has the ad¬ 
vantage of computational simplicity and speed, however 
the resulting layers are not necessarily representative of 
the true physical situation and therefore benchmarking 
against other methods and experiment is critical. Gen¬ 
erating stoichiometry or density gradients across an ar¬ 
tificial junction is not something that can be simulated 
directly using this process, so to investigate the effect of 
these properties, a number of constant density and stoi¬ 
chiometry models were produced. A more sophisticated 
method, closely mimicking the oxygen deposition process 
and examining the effects of layer thickness will be con¬ 
sidered in future work. 


III. MODEL CONSTRUCTION 



FIG. 1. Model of a Josephson junction comprised of alu¬ 
minium (grey) and oxygen (red). Two superconducting re¬ 
gions composed only of aluminium, separated by an amor¬ 
phous AIO1.25 barrier with a density 0.75 times that of corun¬ 
dum. 

To obtain realistic, high precision atomic positions, 
computational models of the junction were created us¬ 
ing a combination of molecular mechanics and Den¬ 
sity Functional Theory (DFT). A 4x4x5 supercell 
of bulk aluminium representing both the top and bot¬ 
tom slabs was relaxed in the DFT code VASF^^iising 
a projector-augmented wave (PAW) potentia l 491501 , ob¬ 
taining a 16.168 x 16.168 x 20.183 A cell. Exchange- 
correlation interactions were evaluated using the PBE 
functional^; a 7x7x7 F centered Monkhorst Pack K 
point mesh and a plane wave cutoff of 250 eV. 

Formation of the amorphous A10 x layers required a 
number of preparation steps to accurately represent ex¬ 
perimental results. The low temperature and pressure 
phase of aluminium oxide (commonly referred to as 
corundum or a-A^Oa) was used as a basis for all the con¬ 
structed junction models. Experimental investigations of 
stoichiometry suggest, in general, an oxygen deficiency 
with oxide O/Al ratios varying between 0.6 and l.#2, 
which are highly dependent on the fabrication process. 


In response to this, we construct models with four sto¬ 
ichiometries: AIOo.8, AIOi.o, AIO 1.25 and AIO 1 . 5 . The 
oxide density may also be an important formation vari¬ 
able. For simplicity we identify oxide density in multiples 
of the (average) corundum density: 4.05 g/cm 3 , and con¬ 
struct junctions with 0.5, 0.625, 0.75, 0.875 & 1.0 density 
multiples for each stoichiometry listed above. A value of 
3.2 g/cm 3 is typicaP^ (which corresponds to a density 
multiple of 0 . 8 ), although theoretical predictions suggest 
altering the density of this barrier may suppress noise 
sources of the junction^. 

Using AIO 1.25 w hh a density multiple of 0.75 as an 
example, a 6 x 6 x 1 supercell of corundum was geometry 
optimised in the software package GULF 5 ^, employing the 
empirical Streitz-Mintmire potential 53 which can capture 
the variable oxygen charge states when present in a pre¬ 
dominantly metallic environment. This capability is par¬ 
ticularly important here, as a Josephson junction has two 
metal-oxide interfaces. This large superstructure was re¬ 
quired due to the trigonal nature of the lattice, as it was 
then cut down such that the xy plane of the bulk alu¬ 
minium slab could be covered. A non-periodic slab of 
corundum measuring 16.168 x 16.168 x 11.982 A was the 
result of this process. Oxygen atoms were randomly re¬ 
moved from the corundum lattice until the appropriate 
stoichiometry of AIO 1.25 was obtained and the cell was 
shortened in the ^-direction to achieve a 0.75 fractional 
multiple of the corundum density. These changes add 
quite a lot of force onto the structure, so a geometry 
optimisation (in GULP) was undertaken at this stage to 
minimise energy contributions. To simulate the oxygen 
deposition phase and generate the amorphous nature of 
these layers, the structure was then annealed using NVT 
molecular dynamics at 3000 K with a 1 fs step size for 
3 ys and quenched to 350 K over a 1.5 ys period. 

The AIO1.25 layer was inserted between two bulk A 1 
supercells described above with 0.5 A of vacuum space on 
each side. The junction was further annealed to simulate 
a metal-metal-oxide interface reconstruction using VASP 
NVT Molecular Dynamics at 300 K until equilibrium was 
reached (approximately 250 ionic steps), then geometry 
optimised using a 2 x 2 x 1 T centered Monkhorst Pack K 
point mesh and a 450 eV plane wave cutoff to obtain the 
final model, depicted in Figure [I] 

For comparison, junctions were also modelled without 
the added computational overhead of DFT by solely em¬ 
ploying GULP and the Streitz-Mintmire potential. The 
construction process of these models matches the proce¬ 
dure above, but interchanges the ab initio optimisations 
of the oxide layer with an empirical framework. 


IV. RESULTS AND DISCUSSION 

To validate our models against experimental observa¬ 
tions, we perform a number of statistical tests to scruti¬ 
nize the structures. First we must ensure that the oxide 
layer of the junctions are in fact amorphous in nature. 
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We employ a projected radial distribution function 


G(r) = lim - ? - 

dr—>o 47 r (N pairs /V) r 2 dr 


(1) 


where r is the distance between a pair of particles, p(r) 
is the average number of atom pairs found at a distance 
between r and r + dr, V is the total volume of the system, 
and iVp a i rs is the number of unique pairs of atoms 54 . This 
function was calculated for each stoichiometry and den¬ 
sity configuration using oxygen as the reference species, 
and aluminium atoms in the amorphous region along 
with the superconducting bulk as the projection species. 
Figure [2] depicts the results of this analysis. 


o 



r (A) 


FIG. 2. Evolution of the oxygen projected radial distribution 
function G(r). Crystalline corundum (thin, gray); Metal— 
metal-oxide interface reconstruction before (dash dotted, 
green) and after (dashed, blue); final optimised geometry 
(thick, red). Inset: Oxygen projected G(r) computed using 
ab initio (VASP) and empirical (GULP) methods, showing no 
statistically significant differences. 


A major peak is visible centered around 1.85 A, which 
corresponds to convolution of the two Al-0 bond dis¬ 
tances, 1.852 A and 1.971 A of the corundum crystal. 
For a crystalline G(r) this peak is deconvolved to two 
delta functions (see Figure [2] and the discussion below), 
where here we see a broadening of the statistics and hence 
differences in neighbour distances: diverging from a crys¬ 
talline form. Moving away from this peak to larger dis¬ 
tance separations, we see the statistics tending toward a 
uniform result similar to what a liquid would produce un¬ 
der this analysis. These two features represent an amor¬ 
phous system quite well, as close range order suggests a 
connection to the crystalline form whilst long range or¬ 
der no longer agrees with such periodic conditions. It’s 
also significant to note that we don’t observe neighbours 
closer than ~ 1.5 A which is a good indication that the 


models do not have non-physical neighbour forces acting 
on atoms. 

Most importantly, this trend is almost uniform across 
all the modeled junctions, which indicate the process out¬ 
lined in Section |III| is capable of producing amorphous 
oxides whilst varying other physical parameters of the 
system. An evolution of the important steps in the pro¬ 
cedure is depicted in Figure [2] 

The corundum G(r) (thin, gray) is a complicated struc¬ 
ture due to the 30 atom unit cell of the crystal, how¬ 
ever it is clear from this figure where much of the amor¬ 
phous structure originates from. Specifically the 1.852 
A and 1.971 A Al -0 bond distance contributions and 
the void in the 2-3 A range. After the melt/quench 
phase of the procedure the lattice still appears liquid- 
like (dash dotted, green). Whilst the quench cycle min¬ 
imises the possibility of atoms position very close to one 
another due to an excess of kinetic energy, it still ap¬ 
pears to exhibit liquid behaviour. This may be a short¬ 
coming of the Streitz-Mintmire potentials ability to cap¬ 
ture the relevant physics, however this is rectified after 
the metal—metal-oxide interface reconstruction is com¬ 
pleted (dashed, blue) using the ab initio methods. Fi¬ 
nally, the geometry optimisation (thick, red) yields a 
smoother 1.85 A peak and recovers some of the void re¬ 
gion around 2 A. 

The inset of Figure[2]compares the optimal G(r) results 
for both the VASP and GULP simulations. Whilst these 
results are very similar, the GULP simulation actually 
produces a drastically different final structure. We find 
under GULP simulation that stoichiometric ratios higher 
than 1:1 are not stable and oxygen atoms diffuse into the 
metallic regions until a stoichiometric ratio of at most 
1:1 is achieved. As a result of this excess oxygen diffu¬ 
sion, the junction width can increase by up to 30% or 
more over the course of the simulation. At high densi¬ 
ties and stoichiometries (higher than typical amorphous 
alumina) some expansion of the oxide region is also seen 
in the ab initio simulations, although this effect is much 
less pronounced. Higher oxygen mobility in GULP could 
be attributed to shortcomings of the empirical potential, 
but we see very little increase in oxide distribution dur¬ 
ing the optimisation phase - suggesting that the details 
of the Nose-Hoover thermostat routine employed during 
the MD simulation may play a role. 

The total energy of a computational model is a good 
indication of the structure’s electronic stability. Due 
to the stoichiometry changes invoked in the oxygen de¬ 
pleted models, not all structures have the same number 
of atoms. This gives structures with more atoms (such 
as AIO 1 . 5 ) additional electronegativity which in turn re¬ 
sults in a deeper potential well and a large total energy. 
In order to be able to validly compare systems of different 
stoichiometry, we normalise the total energy of each sys¬ 
tem by a factor \F\. F = is calculated as the 

linear combination of the number of atoms of chemical 
species k (N k ) by the chemical potential of that species 
(Hk), where k = {Al, O}. The chemical potential for alu- 
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minium, fj ,ai was obtained by calculating the DFT total 
energy of a 4x4x5 supercell of bulk Al and dividing by the 
number of atoms in the supercell. Similarly the chemi¬ 
cal potential of oxygen was obtained from calculating the 
DFT total energy of a 2x2x2 supercell of bulk AI 2 O 3 
using /io = (hai 2 0 3 ~ 2/iAi)/3, where p.Ai 2 0 3 is the total 
energy of a molecular unit of AI 2 O 3 . The factor F (essen¬ 
tially the free energy at T = 0) effectively allows one to 
partition the total energy of the system using the chem¬ 
ical potentials of each component species, as a means to 
compare the energies of systems with differing number of 
chemical components. 

It is clear from Figure [3] that stoichiometry plays a 
larger role in energy minimisation than density, and that 
the structures would prefer additional oxygen to min¬ 
imise internal forces. This suggests that fabrication pro¬ 
cesses that generate oxygen deficiencies may be inviting 
the inclusion of alien species or oxygenic site hopping in 
an attempt to rectify this offset. 

Density changes seem to alter the energy contribu¬ 
tion marginally. Minimum energies correspond to den¬ 
sity multiples between 0.6 and 0.75, slightly lower than 
typical constructions of 3.2 g/cm^H (an 0.8 density mul¬ 
tiple); which may indicate another method of experimen¬ 
tally optimising the junction formation process. 

Coordination number is a useful metric which allows 
for some insight into both the crystallinity of the struc¬ 
tures being analysed, and their similarity to fabricated 
junctions. For instance, in the corundum structure ev¬ 
ery aluminium ion is coordinated with six oxygen ions. In 
amorphous alumina, the proportion of 6 -coordinated alu¬ 
minium as compared to 4-coordinated aluminium is an 
experimentally accessible quantity and has been reported 
on previously®. However, in order to establish this ra¬ 
tio it is assumed that there is a bimodal distribution of 
octahedral (AlCD) and tetrahedral (AIO4) coordination. 
Ratios of A 106 :A 1 C >4 are quoted in a range from 80:20 
to 30:70, depending on the method by which the oxide 
layer was formed 57 . More modern techniques using nu¬ 
clear magnetic resonance (NMR) are also able to resolve 
the AIO5 coordination®. 

Figure [4] shows the distribution of oxygen coordina¬ 
tion about aluminium as a function of density and sto¬ 
ichiometry. These results are calculated using an ADO 
bond length cutoff of 2.5 A, which corresponds to the 
first minimum after the nearest neighbour peak in the 
G(r) (see Figure [ 2 ]). As one would expect, the coordina¬ 
tion number (for ADO bonding) increases with increas¬ 
ing density or stoichiometry. We also note that there 
exists a reasonable proportion of 2- and 3-coordinated 
aluminium atoms, which persists at high density and sto¬ 


ichiometry. In order to compare directly to previous ex¬ 
perimental and theoretical work, we compute the ratio 
of 4-, 5- and 6 -coordination for Al O bonding, matching 
the stoichiometry of 1.5 and assuming the density mul¬ 
tiple closest to experimental values (0.750). The results 
are presented in Table[I] We observe excellent agreement, 
both before and after the ab initio optimisation. 

TABLE I. Relative proportions of 4-, 5- and 6-coordinated 
aluminium atoms within the oxide layer for a density of 0.75 
and stoichiometry of 1.5. 



4(%) 

5(%) 

6 (%) 

VASP (before optimisation) 

57 

39 

4 

VASP (after optimisation) 

53 

43 

4 

Lee et experiment 

55 ±3 

42 ±3 

3 ± 2 

Momida et alW\ theory 

60.4 

29.2 

10.4 


V. CONCLUSIONS 

Precise computational models of Josephson junctions 
are becoming crucial to efforts to reduce dissipation and 
loss in superconducting circuits. The limits of compu¬ 
tational resources mean that full ab initio models are 
computationally intractable. However, a combination of 
ab initio and empirical models holds promise for develop¬ 
ing flexible and efficient simulation approaches. We have 
constructed models of amorphous aluminium-oxide barri¬ 
ers, sandwiched between crystalline aluminium. Through 
comparisons with both previous theoretical analysis and 
experimental measurements, we have shown that the re¬ 
sulting structures are representative of those fabricated 
experimentally. The structure of such junctions can be 
used as input conditions to potential microscopic models 
for either charge or magnetic defects in Josephson junc¬ 
tions. Through this approach, free parameters in existing 
phenomenological defect models can be determined via 
information directly obtained from the atomic positions. 


ACKNOWLEDGEMENTS 

This research was supported under the Australian 
Research Council’s Discovery Projects funding scheme 
(project number DP140100375). Computational re¬ 
sources were provided at the NCI National Facility sys¬ 
tems at the Australian National University through the 
National Computational Merit Allocation Scheme sup¬ 
ported by the Australian Government. 


* tim.dubois@tcqp.science 

1 J. Clarke, A. N. Cleland, M. H. Devoret, D. Esteve, and 
J. M. Martinis, Science 239 , 992 (1988) 


2 J. Martinis, S. Nam, J. Aumentado, and C. Urbina, Phys¬ 
ical Review Letters 89 , 117901 (2002) 

13 B. L. T. Plourde, T. L. Robertson, P. a. Reichardt, 
T. Hime, S. Linzen, C. E. Wu, and J. Clarke, Physical 












5 



FIG. 3. Normalised total energy for various junction models with stoichiometry (left) and density (right). Although the energy 
is strongly dependent on stoichiometry, we see a trend to optimal densities of approximately 75% of the density of corundum. 
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FIG. 4. Distribution of oxygen coordination about aluminium as a function of density and stoichiometry, showing a tendency 
to higher coordination number with increasing density or stoichiometry. 
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